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Abstract. 



^/^ ' We present the complete ground state phase diagram of the Holstein model in two and three 

^ ' dimension considering the phonon variables to be classical. We first establish the overall 

f^ , structure of the phase diagram by using exact diagonalisation based Monte Carlo (ED-MC) 

OO , on small lattices and then use a new "travelling cluster" approximation (TCA) for annealing 

>--^ ' the phonon degrees of freedom on large lattices. The phases that emerge include a Fermi liquid 

>0 ' (FL), with no lattice distortions, an insulating polaron liquid (PL) at strong coupling, and a 

^^~^ I charge ordered insulating (COI) phase around half-filling. The COI phase is separated from the 

\l , Fermi liquid by a regime of phase coexistence whose width grows with increasing electron-phonon 

^p ' coupling. We provide results on the electronic density of states, the COI order parameter, and 

■4.^ ' the spatial organisation of polaronic states, for arbitrary density and electron-phonon coupling. 

Cu ' The results highlight the crucial role of spatial correlations in this strong coupling problem. 
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Introduction. - Electron-phonon (EP) interactions are ubiquitious in metals and dominate 
the finite temperature resistivity in most electron systems. While the perturbative, Fermi 
liquid, regime in EP systems is well understood [1], the electronic ground state is fundamentally 
reorganised when EP interactions are strong. The presence of a strong local interaction can 
C^ I generate a large lattice distortion, creating a potential well for the electron, and lead to a 'self 

trapped' polaronic state [2]. Strong electron-phonon coupling, sometimes in conjunction with 
other interactions, form a crucial component in the physics of several correlated systems, e.g, 
the manganites [3] , the nickelates [4] , and the traditional charge density wave systems [5] . 

The full 'polaron problem', considering the quantum dynamics of the phonons and strong EP 
coupling, is well understood for a single (or few) electrons [6, 7, 8, 9, 10]. The non perturbative 
nature of the problem and the exponential growth in Hilbert space with electron number, Nei, 
has made the finite density problem difficult to access. Usually, the physical interest is in the 
"adiabatic" regime since typical phonon frequencies in most materials are much smaller than 
the electron hopping scale. Although adiabaticity by itself does not lead to a simpler problem. 
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the adiabatic limit, where phonon dynamics is ignored, leads to a relatively tractable situation. 
This limit, of electrons coupled to classical phonons, has been explored within dynamical mean 
field theory (DMFT) [11, 12, 13]. DMFT, however, loses out on spatial correlations or the 
possibility of accessing non periodic phases, important at strong coupling. In this paper we 
use an approach that explicitly retains spatial correlations and set out the ground state of 
many electron systems, in two and three dimension, coupled to adiabatic phonons. 
Let us define the Holstein model, whose adiabatic limit we explore: 

H = -iE^S - mE"« - AE"^^;, + ^Y.p' + f E^? (1) 

(ij) i i i i 

The t are nearest neighbour hopping on a d dimensional lattice, /i is the chemical potential, 
and rii = c\ci is the electron density operator for spinless fermions. The phonon coordinate is 
Xi, and Pi is the momentum conjugate to x^, with [xi,pj\ = iSij. M and K are respectively 
the 'mass' and stiffness of the phononic oscillators, and A is the EP coupling. We set t — I, 
as our reference scale, put K ~ 1 and also h = 1. The adiabatic limit sets M ^ oo. For 
"quantum" phonons, with ujph = y^K/M, the last three terms in H can be written as: 

^Aq J2i "-i(^l +bi) + ^ph J2i ^l^j' "^sing Xi = (bi + b})/ y^2Mujph, and Ag = \y/ujph/i2K). 

The parameter space of the model is defined by: (i) the 'coupling parameter' a = Ep/2dt, 
with Ep = }? jlK, (ii) the 'adiabaticity' 7 = LOph/t, and [iii) electron density, n, in addition 
to dimensionality, d. The parameter a quantifies the competetion between trapping and 
delocalisation, while 7 measures the 'slowness' of the phonons compared to electron dynamics. 

The ground state of this model is understood only in a few limiting cases: (i) At weak to 
moderate coupling, traditional Migdal-Eliashberg (ME) theory [14] can be used for the many 
electron problem for 7 <C 1, and describes a Fermi liquid (FL) metal (or a superconductor, 
when electron spin is also considered), iii) At strong coupling, and for arbitrary 7, there is 
no controlled analytic theory even for a single electron coupled to the phonon system. There 
are, however, powerful numerical methods and the 'single polaron' ground state is essentially 
understood [6, 7, 8, 9, 10]. Little is known of the many electron problem beyond ME theory. 
{Hi) Controlled theory of strong coupling finite density systems has been possible only in the 
classical limit, 7 = 0, in the limit d —^ 00. The original study by Millis ct al. [11] highlighted 
the FL and polaronic insulator phases, within DMFT, and other DMFT studies have explored 
charge ordering [12, 13]. 

While contributing to the initial understanding, DMFT misses out on some crucial aspects 
of the physics, (i) The spatial correlation of lattice distortions, and the resulting correlations 
between polarons, when they form, is not accessible, and (ii) the possibility of phase coex- 
istence and cluster formation is hard to access. We need to move beyond DMFT to recover 
these physical effects, and address current materials issues. The method we use, described 
next, can handle many electron systems strongly coupled to adiabatic phonons, at arbitrary 
disorder, explicitly in 1 — 3 dimension. This paper focuses on the ground state of 'clean' 
systems, subsequent papers discuss the finite temperature effects and the role of disorder [15]. 



Method. - At low temperature the adiabatic EP problem reduces to finding the phonon 
configuration that would minimise the total energy: £{x} = {^t'Yliiij)^i^-j ^ IJ-^iiT-i ^ 
AX^i ''T-iXi){x} + {K/2) ^^ a;?. At strong coupling the electronic energy (in angular brackets) 
in a phonon background, {x}, is not analytically known, since the problem is equivalent 
to that of electrons in an arbitrary strongly fluctuating landscape. The only unbiased way 
to solve the problem is to generate Monte Carlo samples of the phonon configurations and 
accept or reject them using the fermion energy obtained by exact diagonalisation [16]. Such 
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Fig. 1. - Phase diagrams at T = 0. FL indicates Fermi liquid (band metal), the chessboard pattern 
represents a charge ordered insulator with {n,n, ..} modulation, and PL is a positionally disordered 
insulating polaron liquid. Shading indicates regimes of phase separation. The phase diagram are 
obtained through a combination of TCA and analytic estimates: (a) result on 2d, TCA with Lc — 4: 
and system size 24^, (fo) result on 3d with Lc — 4 and system size 8^. The coexistence width at weak 
coupling is estimated as explained in the text. 



exact diagonalisation based Monte Carlo (ED-MC) is feasible only for electronic Hilbert space 
dimension Hpf upto ^ 100 since the computation time t^v grows as Hj^. Although the finite 
size effects are quite strong, the small sizes do provide a feel for the phase diagram. Our 
'reference results' in 2d and 3d are based on ED-MC on 8^ and 4"^ systems. 

To overcome the severe size limitation of ED-MC we have used a 'travelling cluster' approx- 
imation [17] (TCA) for estimating the energy change due to a phonon move. Thus, instead of 
estimating the energy change of the update Xi — > x^, at site R^, as AS = £{x'} — £{x}, by 
diagonalisation of the full Hamiltonian, we obtain an estimate from a 'cluster' Hamiltonian 
constructed by considering L^ (or L^) sites around the reference site R^. The energy of the 
system as a whole of course cannot be approximated by the energy, £ci{Ri}, of the cluster, but 
the change in the energy of the system when a phonon update is made is accurately captured 
by the smaller subsystem. We have broadly benchmarked TCA against ED-MC in an earlier 
paper [17], and also provide several comparisons between the exact and TCA results in this 
paper. Our large size TCA results are obtained by using 4^ clusters on 24^ systems in two 
dimension, and 4'^ clusters on 8^ in three dimension. 

To allow effective annealing we study the system at low finite temperature, T ~ 0.02, with 
varying fx. At equilibrium, at a given /i, we calculate the following: (z) the electron density, 
n(/x), (ii) the structure factor S{q) = (l/iV^) X;y((n,)(nj))e*i-(^'-^^), including 5(7r,7r, ..), 
the order parameter for commensurate charge ordering {{rii) is the quantum average of rii 
in a MC configuration and the outer angular brackets indicate average over configurations), 
(iii) the electronic density of states (DOS), N{ui), (iv) the distribution of lattice distortions, 
as well as MC snapshots of the electron density. Fig.l sets out the ground state phase diagram 
in 2d and 3d, while the detailed indicators, like n{fi), S'(7r,7r, ..), and N{uj), from which the 
phase diagram is constructed, are shown in Figs. 2-3. To access the "weak coupling" part of 
the phase diagram, A < 1.5 say, we have used an analytic method described later. 



Phases. — The two panels in Fig.l show the asymptotic large L phase diagram obtained 
through a combination of TCA and analytic estimates. Our simulations reveal that there 
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Fig. 2. - Panels (a) and (6) show the variation in carrier density with fi computed with ED-MC 
(symbols) and TCA (lines). Panel (a) is for 2d, with ED on 8^ and TCA with 4^ on 8^, while (6) is for 
3d, with ED on 4^ and TCA with 3^ on 4^. Panels (c) and (d) show the COI 'order parameter' with 
varying electron density and coupling. Panel (c) is for 2d, panel (d) is for 3d. Symbols for ED, firm 
lines for TCA, same system size as in (a) — (6). The insets to the lower panels show the COI order 
parameter at n = 0.5. TCA on large sizes does not significantly change the T — > order parameter. 
Panel (a) and (c) share a common legend, as do panel (fe) and (d). 



are primarily three phases in the adiabatic limit, for T ^ 0, in both 2d and 3d. These are 
(i) the FL phase, i.e, the tight binding electron system, without any lattice distortions, (m) an 
insulating polaron liquid (PL) with 'liquid like' short range positional correlations, and a gap 
in the DOS, and (in) charge ordered insulating (COI) phases, with a gap in the DOS, and 
a peak at q = {7r,7r, ..} in S{q). The COI is either a simple "chessboard" phase at n = 0.5, 
or a phase with "defects" off n = 0.5 (and strong coupling). There is also the possibility of 
other charge ordered phases, apart from {tt, tt, ..}, in the vicinity of n = 1/4, 1/8, etc, at strong 
coupling, but their energy difference with respect to the PL phase is very small, making it 
difficult to access them in a low but finite T simulation. There are no "metallic" phases with 
lattice distortions, in contrast to the DMFT results in the adiabatic limit [11]. We discuss the 
phase diagram in terms of the two broad regimes: (i) strong coupling polaronic phases, and 
(a) weak coupling: the COI instability and phase separation. 

(z) Strong Coupling: A naive balance of the 'polaron binding energy' Ep = X^ /2K (for a 
site localised electron), and the lower edge of the tight binding band, —Ef, — 2dt, suggests 
that polaron formation would occur at Ac « ViKdt, implying A^^ ^ 2.82 and Aj?'' ^ 3.46. 
More accurate earlier estimates [7], and our simulations, indicate X^"^ ^ 2.6 and A^^^ ^ 3.3. 
The lowering arises because the polaron is not quite site localised, and taking into account 
the kinetic energy ~ dt'^ /Ep, the polaronic state becomes favourable slightly before the naive 
threshold. The polaron states which emerge for A > Ac are compact, although not site localised, 
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Fig. 3. - Density of states computed with TCA on 24 x 24 in 2d, and 8 x 8 x 8 in 3d. Symbols: circles 
for n = 0.1, triangles for n — 0.3, and squares for n — 0.5. (a) 2d, A — 2.2, (b) 2d, A = 3.0, (c) 3d, 
A = 2.6, (d) 3d, A = 3.5. 



with wavefunctions decaying over a couple of lattice spacings. Due to the strong localisation, 
and the associated short range lattice distortion field, the finite density system is a short range 
correlated array of trapped electrons (the PL phase), keeping maximum mutual separation in 
order to maximise kinetic energy gain by virtual hopping. In the dilute limit the polaron system 
has a very large compressibility, dn/dfx, since the polarons are effectively non interacting and 
H is virtually 'pinned' to the single polaron energy as the carrier density is increased from zero. 
This dictates the sharp rise in the n — fj, curves. Fig. 2. (a) — (5). However, with n increasing 
towards half-filling the short range 'repulsion', ^ 0{t^ /Ep) for nearest neighbours, begins to 
be felt and the n — ii characteristic is flat at n = 0.5, indicative of an incompressible half filled 
state. The polaronic repulsion favours a "chessboard" CO phase at n = 0.5. 

The off half-filling phase maintains the positional correlations of the n — 0.5 phase, but has 
vacancies (or defects), see e.g, Fig. 4. The polarons tend to avoid nearest neighbour locations, 
preferring to be along the diagonal (where the 'repulsion' is weaker). As long as the system is 
sufficiently dense, these positional correlations 'percolate' sustaining long range order, with a 
clear peak in S{tt, tt, ..). Our results suggest that the strong coupling off half-filling COI phase 
survives down to a density « 0.35 in the coupling regime shown (Fig.l). The order parameter 
for the COI phases is shown in Fig.2.(c) — (d), and the DOS in Fig. 3. The order parameter 
computed with TCA on small sizes, see Fig. 2 caption, accurately matches ED-MC results. 
TCA results on large sizes for ^(TrjTr, ..) (not shown) are not significantly different. 

The critical coupling for the FL to PL transition is density dependent, decreasing signif- 



EUROPHYSICS LETTERS 




Fig. 4. - Monte Carlo snapshots of the electron density rir with varying n and A in 2d. Top panels 
A = 2.2, bottom panels A = 3.0. Density n = 0.10, 0.30, 0.50, left to right. Top left panel is FL. Both 
panels to the extreme right are 'perfect' COI, while the other three panels are PL. 



icantly, Fig.l (c) — (d), from the single polaron threshold as n is increased. By the time 
n ~ 0.3 — 0.4, Ac is reduced to « 75% of the 'single polaron' value, in both 2d and 3d, i.e 
the critical Ep is almost halved. Qualitatively, electrons which are already localised affect an 
added particle via their strongly inhomogeneous distortion field, {x}, which interplays with 
EP coupling leading self consistently to the decreasing Xc{n). Trapping at finite n involves a 
combination of (self generated) "disorder" and polaronic tendency [18], akin to the interplay 
of Anderson localisation and single polaron formation. All the phases, COI or PL, which occur 
for A > Xc{n) have a clear spectral gap, and are distinguished mainly by 5(q). 

(a) Weak Coupling: At weak coupling there is an instability in a bipartite tight binding 
lattice, at n = 0.5, due to Fermi surface nesting. The density response function x(q, cj = 0) 
diverges for q -^ {7'',7r, ..}. For A ^ 0, where the analysis in terms of the non interacting 
susceptibility is valid, this leads to a charge order instability, with ordering temperature 
Tco ~ te^*/'^''. With increasing A the CO phase continues to be the ground state, but 
as the charge modulation increases, and the band splitting grows, the FL to COI transition 
at r = becomes first order, with respect n, since the electron density in the modulated Xi 
background is significantly different from that in the homogeneous background. The n — fi 
data in Fig. 2. (a) — (&) highlight the discontinuity in n near half-filling. It is difficult to access 
the instability and the coexistence width at weak coupling in a small system since x(q, a; = 0) 
has strong size dependence. Direct simulation can locate the CO phase only down to A ~^ 1.5. 
At lower A we variationally computed the COI order parameter at n = 0.5 on large lattices, 
and also the /Uc(A) at which the FL to COI transition occurs. It suggests that the FL to COI 
transition becomes (noticeably) first order, with a density discontinuity, for A > 1.0, in both 
2d and 3d. The weak coupling COI phase at n = 0.5 connects continuously to the strong 
coupling "polaron ordered" phase [12, 13]. 



Discussion. - Despite accurate handling of strong coupling and spatial correlations, a few 
physical effects in the adiabatic limit are still difficult to capture within a numerical simulation. 
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Other effects, of possible relevance to real materials, require us to go beyond the adiabatic 
limit itself. This section briefly discusses these issues, (i) As we have discussed, even in the 
adiabatic limit, there are possibly additional commensurate CO phases at n = 1/4, 1/8 etc, 
which could show up at very low temperature. There could also be incommensurate CO phases 
at weak coupling, these are hard to access on small systems, given the strong finite size effects 
in x(q, a; = 0). {ii) In real materials ujph would be finite, and quantum fiuctuations would 
modify some of our conclusions: (a) The simple 'band metal' (FL) obtained at LOph — would 
become a correlated Fermi liquid, with effective mass renornialisation, band narrowing, etc, 
if cOph 7^ 0. (6) Quantum fluctuations would restore translation invariance in the PL phase 
[19] so that below a coherence scale, T^oh ^ '-^pin the resistivity falls quickly to zero, instead 
of diverging as the gapped DOS here would suggest, (c) The COI phase would survive for 
n ^ 0.5 but possibly with reduced ordering temperature. Overall, even with finite cuph, as long 
as LOph/t <C 1, our "ground state" phase diagram would be qualitatively useful for T > Tcoh in 
the model. 

In summary, we have used a new real space method that naturally incorporates the spatial 
correlations vital at metallic densities, to clarify the phase diagram of the adiabatic Holstein 
model. This approach allows unconstrained optimisation, not restricted to any specific kind 
of order, and readily reveals the regimes of phase coexistence. It also extends naturally to 
incorporate the effects of disorder and thermal fluctuations [15]. Another possibility, we will 
separately report, is to include spin degrees of freedom, via double exchange, to model the 
phenomena in manganite physics. 
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